* ********************************************************************************************
* The Effects of Weather Shocks on Economic Activity: What are the Channels of Impact?
* Sebastian Acevedo, Mico Mrkaic, Natalija Novta, Evgenia Pugacheva, Petia Topalova
* With support from Gavin Asdorian, Olivia Ma, Jilun Xing and Yuan Zeng 
* Replication files for Table 3; Figure 7; panel 1-2
* ********************************************************************************************
clear all
set more off
set matsize 11000

cd "C:\Users\..\replication" // set your working directory to where the replication folder is saved

* the Y variable that will be used in regressions
local Y_list ava_rlcu_wdi ag_prd_crop mva_rlcu_wdi sva_rlcu_wdi ni_r nm_r mrt_i hdi

* Horizons
global k 7
* ****************************************************************************************
* Data
* ****************************************************************************************
use "data/country_dataset.dta", clear
keep ifscode country year GGDC_* wdi_region ae em lics tmp_UDEL* tmp_CRU* precip_UDEL* precip_CRU* 

keep if GGDC_AGR_VA!=.		

foreach X in AGR MIN MAN PU CON WRT TRA FIRE GOV OTH {
	gen GGDC_p_`X' = GGDC_`X'_VA_Q05/GGDC_`X'_EMP
	cap rename GGDC_`X'_EMP emp`X'
	cap rename GGDC_`X'_VA va`X'
	cap rename GGDC_p_`X' prod`X'
	}


reshape long emp va prod, i(ifscode country year wdi_region tmp_CRU_pw50 precip_CRU_pw50 em ae lics) j(sector) string


* create variable for region-sector-year, encode it to be numeric, to be used as fixed effects
gen tempvar1 = wdi_region + " " + sector+" " + string(year) if wdi_region != ""
encode tempvar1, generate(ry)
drop tempvar1

* generate squared terms
foreach X in UDEL_pw10 UDEL_pw90 UDEL_pw50 CRU_pw90 CRU_pw10  CRU_pw50 {
	gen tmp_`X'2 = tmp_`X' ^ 2
	gen precip_`X'2 = precip_`X' ^ 2
}

gen id_prod = string(ifscode)+" "+sector
encode id_prod, gen(idprod)

tsset idprod year

* calculate Y log difference at different horizons
foreach Y of varlist emp prod va {
	gen ln_`Y' = ln(`Y')
	bys idprod (year), sort: gen d0_ln_`Y' = ln_`Y' - l.ln_`Y'
	forval t = 1/7 {
		bys idprod (year), sort: gen d`t'_ln_`Y' = f`t'.ln_`Y' - l.ln_`Y'
	}
}

**** OUTDOOR SECTORS ****
gen out= (sector=="AGR" | sector=="CON" | sector=="MAN" | sector=="MIN" | sector=="TRA" | sector=="PU")
foreach X in tmp_CRU_pw50 tmp_CRU_pw502 precip_CRU_pw50 precip_CRU_pw502 {
	gen out_`X' = out*`X'
	gen ind_`X' = (1-out)*`X'
	}

	
xi i.idprod	i.ry

* ****************************************************************************************
* Program
* ****************************************************************************************

gen period = -1 in 1
foreach X in ind out {
	gen t_`X'_coef = 0 in 1
	gen t_`X'_se = 0 in 1
	gen t_`X'_pvalue = 0 in 1
	gen t2_`X'_coef = 0 in 1
	gen t2_`X'_se = 0 in 1
	gen t2_`X'_pvalue = 0 in 1
	
	gen p_`X'_coef = 0 in 1
	gen p_`X'_se = 0 in 1
	gen p_`X'_pvalue = 0 in 1
	gen p2_`X'_coef = 0 in 1
	gen p2_`X'_se = 0 in 1
	gen p2_`X'_pvalue = 0 in 1
	
	foreach G in ae em lic {
		* temperature
		gen t_`X'_`G'_coef = 0 in 1
		gen t_`X'_`G'_se = 0 in 1
		gen t_`X'_`G'_pvalue = 0 in 1
		* precipitation
		gen p_`X'_`G'_coef = 0 in 1
		gen p_`X'_`G'_se = 0 in 1
		gen p_`X'_`G'_pvalue = 0 in 1
	}
}
gen nobs = 0 in 1
gen ncountry = 0 in 1
gen r2 = 0 in 1

cap program drop lpmregsctr
program define lpmregsctr
syntax,  yy(string) tt(string) pp(string)
forval t = 0/$k {
	areg d`t'_ln_`yy' out_`tt' out_`tt'2 out_`pp' out_`pp'2 l.out_`tt' l.out_`tt'2 l.out_`pp' l.out_`pp'2 ///
				 ind_`tt' ind_`tt'2 ind_`pp' ind_`pp'2 l.ind_`tt' l.ind_`tt'2 l.ind_`pp' l.ind_`pp'2 ///
				l.d0_ln_`yy' _Iidprod*, cluster(ifscode) absorb(ry)


	replace period = `t' in `=`t'+2'
	
	foreach X in ind out {
		*** Temperature
		replace t_`X'_coef = _b[`X'_`tt'] in `=`t'+2'
		replace t_`X'_se = _se[`X'_`tt'] in `=`t'+2'
		replace t_`X'_pvalue = 2 * ttail(e(df_r),abs(_b[`X'_`tt']/_se[`X'_`tt'])) in `=`t'+2'
		replace t2_`X'_coef = _b[`X'_`tt'2] in `=`t'+2'
		replace t2_`X'_se = _se[`X'_`tt'2] in `=`t'+2'
		replace t2_`X'_pvalue = 2 * ttail(e(df_r),abs(_b[`X'_`tt'2]/_se[`X'_`tt'2])) in `=`t'+2'
		*** Precipitation
		replace p_`X'_coef = _b[`X'_`pp'] in `=`t'+2'
		replace p_`X'_se = _se[`X'_`pp'] in `=`t'+2'
		replace p_`X'_pvalue = 2 * ttail(e(df_r),abs(_b[`X'_`pp']/_se[`X'_`pp'])) in `=`t'+2'
		replace p2_`X'_coef = _b[`X'_`pp'2] in `=`t'+2'
		replace p2_`X'_se = _se[`X'_`pp'2] in `=`t'+2'
		replace p2_`X'_pvalue = 2 * ttail(e(df_r),abs(_b[`X'_`pp'2]/_se[`X'_`pp'2])) in `=`t'+2'
	}
	
	replace nobs = `e(N)' in `=`t'+2'
	replace ncountry = `e(N_clust)' in `=`t'+2'
	replace r2 = `e(r2_a)' in `=`t'+2'
	
	foreach X in ind out {
		*** Temperature
		lincom _b[`X'_`tt']+2*11*_b[`X'_`tt'2]
		replace t_`X'_ae_coef = r(estimate) in `=`t'+2'
		replace t_`X'_ae_se = r(se) in `=`t'+2'
		replace t_`X'_ae_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'
		
		lincom _b[`X'_`tt']+2*22*_b[`X'_`tt'2]
		replace t_`X'_em_coef = r(estimate) in `=`t'+2'
		replace t_`X'_em_se = r(se) in `=`t'+2'
		replace t_`X'_em_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'
		
		lincom _b[`X'_`tt']+2*25*_b[`X'_`tt'2]
		replace t_`X'_lic_coef = r(estimate) in `=`t'+2'
		replace t_`X'_lic_se = r(se) in `=`t'+2'
		replace t_`X'_lic_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'
		
		*** Precipitation
		lincom _b[`X'_`pp']+2*8*_b[`X'_`pp'2]
		replace p_`X'_ae_coef = r(estimate) in `=`t'+2'
		replace p_`X'_ae_se = r(se) in `=`t'+2'
		replace p_`X'_ae_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'
		
		lincom _b[`X'_`pp']+2*9*_b[`X'_`pp'2]
		replace p_`X'_em_coef = r(estimate) in `=`t'+2'
		replace p_`X'_em_se = r(se) in `=`t'+2'
		replace p_`X'_em_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'
		
		lincom _b[`X'_`pp']+2*11*_b[`X'_`pp'2]
		replace p_`X'_lic_coef = r(estimate) in `=`t'+2'
		replace p_`X'_lic_se = r(se) in `=`t'+2'
		replace p_`X'_lic_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'
	}
	}
end		


* ****************************************************************************************
* Run (transportation and utilities is outdoor)
* ****************************************************************************************
lpmregsctr, yy(prod) tt(tmp_CRU_pw50) pp(precip_CRU_pw50)
export excel period - r2 in 1/9 using "output/Table_2-3.xlsx", sheet("GGDC") sheetreplace firstrow(var)


